Accurate and realistic initial data for black hole-neutron star binaries. 
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This paper is devoted to the computation of compact binaries composed of one black hole and one 
neutron star. The objects are assumed to be on exact circular orbits. Standard 3+1 decomposition 
of Einstein equations is performed and the conformal flatness approximation is used. The obtained 
system of elliptic equations is solved by means of multi-domain spectral methods. Results are 
compared with previous work both in the high mass ratio limit and for one neutron star with very 
low compactness parameter. The accuracy of the present code is shown to be greater than with 
previous codes. Moreover, for the first time, some sequences containing one neutron star of realistic 
compactness are presented and discussed. 
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I. INTRODUCTION 



Motivated by the various gravitational wave detectors coming online [l|, 0| , numerical simulations of binary compact 
objects have been extensively investigated in the last years. With progress made for the evolution of both binary 
' neutron stars (BNS) [|| and binary black holes (BBH) @, H, Q , it seems timely to turn to the last interesting type of 
binary: a system composed of one black hole and one neutron star (BHNS). The evolution synthesis codes indicate 
t^J- ' that the detection rate of BHNS with LIGO/ Virgo will be as high (if not higher) than for BNS systems (see Table 
6 of 7]). Simulations of coalescing binary systems usually proceed in two steps. First, one needs to produce initial 
data that verify the Einstein constraint equations and that are as physically relevant as possible. Then, the initial 
configurations are evolved forward in time. Those steps involve rather different techniques and are equally challenging. 

Most of the initial data for coalescing binaries rely on the quasiequilibrium hypothesis that assumes that the objects 
are on exact closed circular orbits. This is of course only an approximation because no closed orbits can exist for 
those systems in general relativity. However, this description should be rather accurate, at least for lar ge s eparations. 
^ { The quasiequilibrium approximation has been applied to BNS [1, [l(| [HI, [12] and BBH systems [HI, 1T3 , [l5l . [lfjj . In 
the BNS case, the evolution of such initial data have exhibited circular- like trajectories [|[. The long term evolution 
£~ i ■ of the data in the BBH case is still underway. The aim of this paper is to apply the same technique to a BHNS, 
without assuming extreme mass ratios like in [l7|]. On the contrary, a moderate mass ratio of 5 is used. This value 
. is chosen mainly for comparisons with previous work [l8j but is sufficient to demonstrate the ability of the code to 
handle binaries with close masses. Let us mention that the simulations of the formation of compact binaries indicate 
a slightly lower mass ratio to be more probable [l9j]. The present work shares some properties with [ijj] even if 
some details are different. Moreover, one of the main shortcomings of [l8| is the fact that only a NS of unrealistic 
compactness, as small as 0.0879, is considered. This implies a mass close to 0.7 M© for most of the available EOS, 
which is much smaller than the canonical value of 1.2 — 1.4Mq. On the contrary, in this work and for the first time, 
a BHNS with a realistic NS is computed. Indeed, the most compact star presented here has a compactness of 0.15 
which makes its mass in the range 1.2 — 1.5 M depending on the EOS. 



II. EQUATIONS 

The standard 3+1 decomposition of the Einstein equations is used, in which the spacetime is foliated by a family 
of spatial hypersurfaces. The 4-dimensional metric is given in terms of the lapse function N, the shift vector j3 l and 
the spatial metric Yy . The evolution of the spatial metric is described by the second fundamental form: the extrinsic 
curvature tensor which is the Lie derivative of 7y along the normal to the hypersurfaces. 

The conformal factor '5 is defined by 7^ = , E r 7^ and by demanding that the determinant of the conformal metric 
jij is 1 (so that the determinant of 7^ is \1/ 12 ). A similar conformal decomposition is also performed on the extrinsic 
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curvature tensor. 

Before solving the constraint equations, there are some quantities that must be chosen. Those are the so-called 
"freely specifiable" variables. By this, one means that a solution of the constraint equations can be found for every 
choice of those variables. In the standard "thin-sandwich" approach, these variables are: the trace of the extrinsic 
curvature K and its time derivative dtK, the conformal metric 7y and its time derivative Uij = dtjij. By working in 
the corotating frame, the quasicquilibrium approximation amounts to neglecting the variations with time of all the 
quantities. In particular, one sets dtK = and iy = 0. Maximum slicing K = and conformal flatness 7^ = /y, /y 
being the 3-dimensional flat metric, are also assumed. Even if the spatial metric for a binary system can not be flat, 
this approximation has proven to be more accurate that one could expect for both BBH and BNS 0, [2l], [22j . Let 
us note that the choices for K and 7^ are different from the ones made in [l8| where the authors used values derived 
from the Kerr-Schild metric. 

The mathematical problem is then to solve the following set of five coupled elliptic equations for the two scalars W 
and N and the vector field (3 l : 

AN = 4tt7V* 4 (E + S) + A - 2D, In ^&N (1) 

A* = -2ttV 5 E AijA^ (2) 

8 

A/3 1 + )-D l D P = 16tt7V* 4 (E + p)W + 2A lJ (D N - §ND 3 In (3) 

where all the operators are associated with the flat metric. The conformal extrinsic curvature tensor = ^K 1 ^ 
relates to the lapse and shift vector by A lJ = 1/2N (D l [3' J + & ft - 2/3D k p k f j ). E, S, p and U l are matter terms 
describing the NS fluid: p is the pressure, and E and S are the matter energy density and the trace of the stress 
energy tensor, both measured by the Eulerian observer (whose 4-velocity coincides with the normal to the spatial 
hypersurfaces). U l is the fluid 3- velocity with respect to the Eulerian observer. Explicit expressions for all those 
matter terms can be found, for example, in [l2T ]. In this paper, a polytropic equation of state relates the pressure to 
the fluid baryon number density by p — nnF . In all the following T is fixed to 2 and k is varied to construct NS with 
various compactness. The star is also assumed to be irrotational. It implies that the velocity U l relates to a potential 
$ via : 

if = — r- — L> 4 $, (4) 

where h is the fluid specific enthalpy and r„ the Lorentz factor between the Eulerian observer and the fluid (see [l2[ 
for explicit expressions). The potential obeys the following equation : 

£i?A$ + D l HD^ = ^ 4 hT n U^D t H + £H (D l ^D t (H - /?) + ^hU^DiTn) , (5) 

where H = hxh, £ = d\nH/dln.n, (3 = In (^ 2 7V) and Uq is the 3- velocity of the corotating observer. The irrotational 
fluid also admits an integral of motion which can be expressed as : 

r 

hN — = const., (6) 
To 

where T is the Lorentz factor between the corotating observer and the fluid and Tq between the corotating and the 
Eulerian observers. 

The black hole is described in the framework of the apparent horizons, where non-trivial boundary conditions are 
imposed on a sphere, (see the work of [23| or [24| for a review). This idea has already been successfully applied 
to BBH systems (see [l6[ for the last application to date). However, in this work, like in [TJ, [l5j], the lapse is set 
to on the horizon. In doing so one needs to make a small correction on the shift vector to get a regular extrinsic 
curvature tensor on the horizon. In the case of the BBH it has been shown that this correction was small enough 
[T5I ]. In the BHNS case, the effect of the neutron star on the horizon is even smaller, making this relative correction 
very small (the relative difference between the original and the corrected shift is at most 2 • 10 -4 for the innermost 
configurations). An important difference with respect to what was done in [l5| concerns the rotation state of the 
black hole. Indeed, it seems unlikely that it will be synchronized with the orbital motion and an irrotational black 
hole will be considered. Like in Sec. VB of [l6| . a local rotation rate is imposed. This is done by imposing that, on 
the horizon, the corotating shift is /3* = fl r (d/dip)\ the angle if being associated with the horizon. In the case of the 
corotating black holes of [15| . f2 r is zero. When considering irrotational black holes, fl r must be determined to ensure 
that the quasi-local spin of the black hole vanishes. The appropriate value of O r must be obtained numerically. This 
is to be contrasted with what is done in [l8[ where the authors impose irrotationality only to the first order (see Sec. 
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VA of [l6|), i.e. by demanding that fl r coincides with the orbital velocity 17q. As will be seen later, the correction to 
this, measured by the ratio f r = fi r /f2o 1S °f the order of 10%. f r is one when the first order is valid, that is when 
mass ratio is infinite, and smaller than unity when the masses become comparable. 

Standard asymptotic flatness is used to get appropriate boundary conditions at infinity. 

III. NUMERICS, TESTS AND COMPARISONS 

The system ([T]G|]) is solved using the LORENE library [25J]. This library is developed mainly at the Meudon site 
of Paris observatory and has been successfully applied to the computation of various problems in general relativity 
(see references at [25|). The basic features of Lorene are the following: spectral methods in spherical coordinates 
are employed, mainly using spherical harmonics or trigonometrical functions for the angles (8, tp) and Chebyshev 
polynomials for the radial coordinate r. Space is decomposed in various spherical-like shells, the spectral expansion 
being done in each of those domains. Space is compactified by means of the variable u = 1/r in the outermost domain. 
This enables one to impose exact boundary conditions at infinity, the computational domain covering the whole space. 
Two sets of such domains are used, one centered around each compact object and the equations ([T][3]) are split on 
those two sets of domains (see for example [HJ and [H[ for explicit implementations of this method). For the NS, 
the first domain is adapted to match the surface of the star, thus getting rid of any Gibbs phenomenon that would 
be caused by discontinuities at the surface. With such methods, solving elliptic equations amounts to inverting some 
matrices. Explicit examples are given in [2(| for Poisson equations and the algorithm used to solve Eq. ([5]) can be 
found in Let us mention that the code of [H| is also based the LORENE library. However the two codes have 
been written completely independently and the details of the two implementations are not the same. 

The system ([T][3]), along with the integral Q and the equation §5§ for the potential is solved by iteration, until 
the fields converge to a given threshold (typically 10~ 7 ). During the course of the iteration, some quantities are 
changed in order to fulfill various requirements. For example, the central enthalpy of the NS is rescaled at each step 
in order to get a neutron star of given baryon mass (see Sec. IVD3 of [HJ ). The same technique is also used to 
modify the radius of the black hole to get a given irreducible mass, the position of the rotation axis to ensure that the 
total linear momentum vanishes and the local rotation rate of the black hole to impose irrotationality. The distance 
between the two objects can also be modified by these means. As in [jjj], the orbital velocity is determined so that 
the zero of the gradient of enthalpy lies at the origin of the grid associated with the NS. 

For both objects, space is decomposed on typically 8 or 9 domains, each of them being covered by N r x N# x N v = 
33 x 21 x 20 points. The achieved precision can be measured by some global checks. For instance, if Eq. ^ is fulfilled, 
the ADM mass can be computed in two ways: either by the standard integral at infinity, or by a volume integral 
plus an integral on the BH horizon (see Sec. HID of [11]). For all our configurations, the relative difference between 
the two values is less than a few times 10~ 5 . Other criteria of this kind can be derived that test the other equations 
(using the total angular momentum J and the generalized Smarr formula ; see Sec. IIIC and HID of [Hj]). However, 
in those cases, the precision is also limited by the regularization of the shift. As already stated, the relative difference 
between the original and the corrected shift is at most a few times 10~ 4 which implies a relative difference of 10~ 2 
between the two expressions of J (analogues of Eqs.(67) and (68) of [IB]). The two orders of magnitude between the 
regularization of the shift and the error on the global quantity J is similar to what is observed in the BBH case (see 
Fig. 10 of [H). 

Another test is provided by comparing the results from this work with the ones obtained in the extreme mass ratio 
in [13]. More precisely, in Fig. [TJ the value of \ is shown as a function of the orbital frequency, for a mass ratio 
Mg^/M^g = 10. Mgn denotes the irreducible mass of the BH and M^ s the baryonic mass of the NS, both quantities 
being kept constant along a sequence, x is a measure of the deformation of the star. It is one for a spherical star and 
zero at the mass-shedding limit (see Eq. (52) of [13] for the precise definition). The first panel of Fig. [T] shows the 
value of x f° r -^ns = 0.05 and the second one for M^s = 0.01, the masses being expressed in standard polytropic 
units. The agreement with [13] is good, especially for the isotropic background. The difference is more important 
with the Kerr-Schild approach. 

Another comparison can be made with the sequence presented in [l8j that relies on an approach similar to the work 
presented here. Let us first mention that the neutron star presented in [H| is not realistic, having a compactness as 
small as 3 = 0.0879. For comparison purposes, a sequence with the same parameters has been computed with the 
code described here. In the first panel of Fig. [2J the value of \ as a function of HMq is shown. In all the following, 
Mq denotes the total gravitational mass for infinite separation, which is the sum of the BH irreducible mass Mgg and 
the gravitational mass (i.e. the ADM mass) of the isolated NS: M^ r g av0 . One of the most striking features of the first 
panel is that the data from [l8| behave in a different manner for large separations. Some of this surely comes from 
the Kerr-Schild approach, but the fact that the curve exhibits a very sharp maximum is puzzling. Moreover, for large 
separations, x does not seem to converge to one as it should. On the contrary, the curve from this work exhibits a 
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FIG. 2: Comparison between this work and [18J for a NS of small compactness (H = 0.0879). The ratio between the BH mass 
and the NS mass is 5. The first panel shows the deformation parameter \ an d the second one the binding energy, both as a 
function of the orbital velocity. 



very smooth behavior and clearly goes to one for small values of the orbital frequency. For large values of £IMq, that 
is, for small separations, the two curves agree reasonably well. The second panel of Fig. [2] shows the binding energy 
of the binary defined as = M \dm /Mq — 1, as a function of SIMq. The results from the two approaches are very 
different, the configurations from [l8| being more bounded by a factor of 3. For comparison, the binding energy from 
both Newtonian and 3.5 PN theories are shown (given by Eq. (194) of 27}). The PN result is clearly much closer to 
the result of this paper. This closeness is a strong indication that the conformal flatness approximation is rather good 
because PN theory does not make us of it. The disagreement with [l8| may come from the differences in the choices 
of the "freely-specifiable" variables. If this is the case, it may indicate that the data generated with the Kerr-Shild 
approach have a spurious gravitational wave content. The total angular momenta coming from the two approaches 
are also compared and the agreement is good, the difference being only of the order of 4% which is also the order of 
the difference with the PN result. 
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FIG. 3: Deformation parameter \ (first panel) and binding energy (second panel) as a function of orbital velocity for four 
different compactness parameters S. M^ r s av0 is the same for all the stars and the mass ratio with respect to the BH irreducible 
mass is 5. 



IV. MORE REALISTIC NEUTRON STARS 



As already stated, the NS presented in the previous section has a very small compactness. In order to get more 
realistic neutron stars, this parameter should be increased. This is easily done by decreasing the parameter k in 
the EOS. The baryon mass of the NS is fixed so that they all have the same gravitational mass when isolated, i.e. 
the same Mj|g av . Doing so, BHNS with four different compactness parameters (0.075, 0.100, 0.125 and 0.150) are 
constructed. As in the previous section, the ratio M^/M^ v0 is set to 5. 

The first panel of Fig. [3] shows the deformation x as a function of flMo. As can be expected, the more compact 
stars are less easily deformed and can survive closer to the black hole without being tidally destroyed. The second 
panel shows the binding energy of the binaries, along with the 3.5 PN result for point masses. Once again, the more 
compact the NS, the closer it can get to the BH. It turns out that for the most compact star, the binding energy attains 
a minimum before the NS is destroyed. This may indicate that the system can reach dynamical instability, even if 
true stability can only be investigated by means of dynamical evolutions. The minimum is also marginally attained 
for a compactness of 0.125. The two less compact stars, however, are destroyed before reaching the configuration of 
minimum binding energy. This explains why all the configurations of [l8T | are found to be stable: the compactness of 
the star is small enough that it is destroyed before reaching the innermost stable orbit. The nature of the end point 
of a sequence thus depends on the compactness of the star, and this is an effect that could have some implications 
on the emitted gravitational signal. Fig. Q] shows the value of the lapse (first panel) and of A xx (second panel) in the 
orbital plane for a NS of compactness 0.15. The configuration is the one corresponding to the minimum of binding 
energy. The surfaces of both the NS (thick line on the right) and the BH (thick line on the left) are shown, and one 
can see that the NS is not yet noticeably deformed. 

The total angular momentum J can also be computed. When it admits a minimum (i.e. for the two most compact 
stars) , its position in terms of frequency is consistent with the one for the binding energy (even if they do not coincide 
exactly). The relative importance of the BH local rotation can also be investigated by computing the ratio f r = f2 r /f2o- 
It turns out that f r is close to one (at most 0.91) which shows that the presence of the neutron star has a moderate 
influence on the structure of the BH horizon. f r is almost independent of the compactness of the star. However its 
dependence with frequency is different from Eq. (58) of [l6j |. This is probably simply an effect of the mass ratio 
(5 in this work and 1 in [ill). Finally the virial theorem states that the Komar-like mass and ADM mass should 
be equal for circular orbits. Contrary to the case of BBH [3, HUGH, this is not enforced exactly, the value of the 
orbital velocity being given by equilibrium of the NS fluid. However, the virial theorem can be used as a test of the 
code. It turns out that the virial theorem is verified to better than 2% for all configurations, the dependency with the 
compactness of the NS being moderate. The difference between the two masses is greater for tighter configurations. 
This behavior is similar to what is found in [l8T ] even if the curves do not match, probably because of the different 
choice of "freely specifiable" variables. The violation of the virial theorem must be a measure on how the computed 
configurations differ from exact circularity. In a sense, it reflects the true nature of the movement: a slow inspiral. 

All the configurations presented in this paper have been made public on the LORENE website [25| . 
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V. ERRATUM 



The version of the code used to compute the quasiequilibrium configurations presented in the previous sections was 
faulted by a wrong sign in some of the terms coupling the neutron star to the black hole. Indeed, terms like DiN are 
split in a black hole and a neutron star part DiN\ BH + 75iiV| Ng . The terms generated by the neutron star had the 
wrong sign for the X and Y components, when computed on the grid associated to the black hole. 

With the correct sign, the precision of the computation, as measured by the virial theorem, has improved by at 
least one order of magnitude. The obtained results are now in much better agreement with both post-Newtonian 
results and independent computations by Taniguchi et al. [HI]. This is clearly seen in Fig. [5] which is the correct 
version of Fig. 3. In particular, the more compact the star, the closer the binding energy is to the post-Newtonian 



curve, as expected. The new configurations can be accessed via the Lorene website 25] 
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